#####Replication codes for: Ege, J,. Bauer, M., Wagner, N. and E. Thomann (2022). When does Bureaucracy Matter in the Making of Global Public Policies? Conditions of Policy Influence of International Public Administrations. Governance, DOI: 10.1111/gove.12741.

# install packages
install.packages(c("QCA", "SetMethods", "psych","dplyr", "ggplot2", "ggthemes", "gridExtra"), dependencies = TRUE)

# load packages
library(QCA); library(SetMethods);library(psych); library(dplyr); library(ggplot2); library(ggthemes); library(gridExtra)

# set working directory

# clean workspace
rm(list = ls())

# load dataset

mydata <- read.csv("data_210701.csv", header = TRUE, row.names = 1)
mydata

##### descriptive stats#####
des <- describe(mydata)
des
stargazer(mydata, type="html", summary=TRUE, out="des.html")

##test skewness
skew.check(mydata)

##There are 11 cases of INFL and 6 cases of ~INFL to explain.

##histogram of data

INFL <- ggplot(data=mydata, aes(INFL)) +
  geom_histogram(binwidth = 0.33, fill="lightblue2",
                 colour="lightsteelblue4") +
  scale_x_continuous(breaks = c(0, 0.33, 0.66, 1.0)) +
  scale_y_continuous(name = "Frequency",
                     limits = c(0,8), breaks = c(0, 1, 2, 3, 4, 5, 6, 7, 8)) +
  theme_classic() +
  labs(x="Set Membership", y="Frequency") +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=11)) +
  ggtitle("INFL") +
  theme(plot.title = element_text(hjust = 0.5, size =15, face = "bold"))

AWILL <-   ggplot(data=mydata, aes(AWILL)) +
  geom_histogram(binwidth = 1, fill="lightblue2",
                 colour="lightsteelblue4") +
  scale_x_continuous(breaks = c(0, 1.0)) +
  scale_y_continuous(name = "Frequency",
                     limits = c(0,12)) +
  theme_classic() +
  labs(x="Set Membership", y="Frequency") +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=11)) +
  ggtitle("AWILL") +
  theme(plot.title = element_text(hjust = 0.5, size =15, face = "bold"))


AACTION <-  ggplot(data=mydata, aes(AACTION)) +
  geom_histogram(binwidth = 1, fill="lightblue2",
                 colour="lightsteelblue4") +
  scale_x_continuous(breaks = c(0, 1.0)) +
  scale_y_continuous(name = "Frequency",
                     limits = c(0,14), breaks = c(0, 2, 4, 6, 8,10,12,14)) +
  theme_classic() +
  labs(x="Set Membership", y="Frequency") +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=11)) +
  ggtitle("AACTION") +
  theme(plot.title = element_text(hjust = 0.5, size =15, face = "bold"))


COMPLEX <- ggplot(data=mydata, aes(COMPLEX)) +
  geom_histogram(binwidth = 0.33, fill="lightblue2",
                 colour="lightsteelblue4") +
  scale_x_continuous(breaks = c(0, 0.33, 0.66, 1.0)) +
  scale_y_continuous(name = "Frequency",
                     limits = c(0,8), breaks = c(0, 1,2,3,4,5,6,7)) +
  theme_classic() +
  labs(x="Set Membership", y="Frequency") +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=11)) +
  ggtitle("COMPLEX") +
  theme(plot.title = element_text(hjust = 0.5, size =15, face = "bold"))


CONTEST <- ggplot(data=mydata, aes(CONTEST)) + 
     geom_histogram(binwidth = 1, fill="lightblue2",
     colour="lightsteelblue4") +
     scale_x_continuous(breaks = c(0, 1.0)) +
     scale_y_continuous(name = "Frequency",
    limits = c(0,12)) +
     theme_classic() +
     labs(x="Set Membership", y="Frequency") +
     theme(axis.text=element_text(size=12),
          axis.title=element_text(size=11)) +
     ggtitle("CONTEST") +
     theme(plot.title = element_text(hjust = 0.5, size =15, face = "bold"))
  


gl <- grid.arrange(INFL, AWILL, AACTION, COMPLEX, CONTEST, ncol=3)


#### analysis of necessity#####
colnames(mydata)

# necessary conditions for INFLUENCE
superSubset(mydata, outcome = "INFL", 
            conditions = "AWILL, AACTION, COMPLEX, CONTEST", 
            incl.cut = 0.9, cov.cut = 0.5, ron.cut=0.5)

# necessary conditions for ~INFLUENCE
superSubset(mydata, outcome = "~INFL", 
            conditions = "AWILL, AACTION, COMPLEX, CONTEST", 
            incl.cut = 0.9, cov.cut = 0.5, ron.cut=0.5)

# no necessary conditions meet the thresholds

# results necessity for single conditions

conds <- subset(mydata, select = c("AWILL", "AACTION", "COMPLEX", "CONTEST"))

nec1 <- pof(conds, "INFL", mydata, relation = "nec")
nec2 <- pof(1-conds, "INFL", mydata, relation = "nec")

nec3 <- pof(conds, "~INFL", mydata, relation = "nec")
nec4 <- pof(1-conds, "~INFL", mydata, relation = "nec")

nec1
nec2
nec3
nec4

# produce tables

write.csv(nec1$incl.cov, "nec1.csv")
write.csv(nec2$incl.cov, "nec2.csv")
write.csv(nec3$incl.cov, "nec3.csv")
write.csv(nec4$incl.cov, "nec4.csv")

# plot ~AWILL for ~INFL

xy.plot("~AWILL", "~INFL", mydata, labs = rownames(mydata), necessity=TRUE)

# although ~AWILL has a high consistency necessity, the set is skewed and therefore it is a trivial necessary condition--although there are no deviant cases consistency in kind

###analysis of sufficiency####

#####INFL#####

# truth table for INFL, with all cases
ttOUTCOME <- truthTable(data=mydata, outcome = "INFL", 	
                        conditions = "AWILL, AACTION, COMPLEX, CONTEST", 
                        incl.cut=0.75, sort.by="incl, n", complete=TRUE, show.cases=TRUE)
ttOUTCOME

# truth table for INFL, with dcc only
ttOUTCOME1 <- truthTable(data=mydata, outcome = "INFL", 	
                        conditions = "AWILL, AACTION, COMPLEX, CONTEST", 
                        incl.cut=0.75, sort.by="incl, n", complete=TRUE, show.cases=TRUE, dcc=TRUE)
ttOUTCOME1

# we set the raw consistency threshold at .830, to include those rows with no dcck
ttOUTCOME <- truthTable(data=mydata, outcome = "INFL", 	
                        conditions = "AWILL, AACTION, COMPLEX, CONTEST", 
                        incl.cut=0.830, sort.by="incl, n", complete=TRUE, show.cases=TRUE)
ttOUTCOME

# export truth table

write.csv(ttOUTCOME$tt, "ttINLF.csv")

#conservative solution

csOUTCOME <- minimize(ttOUTCOME, details=TRUE, show.cases=TRUE, use.tilde=TRUE)
csOUTCOME

# produce table
stargazerSol(results = csOUTCOME, outcome = "INFL",
             type = "html", show.cases=TRUE, out = "csOUTCOME.html")


# parsimonious solution
psOUTCOME <- minimize(ttOUTCOME, include="?", details=TRUE, show.cases=TRUE )
psOUTCOME

# produce table
stargazerSol(results = psOUTCOME, outcome = "INFL",
             type = "html", show.cases=TRUE, out = "psOUTCOME.html")


# intermediate solution. We expect that ceteris paribus, AWILL, AACTION, COMPLEX would lead to INFLUENCE. There is no expectation for CONTEST.

isOUTCOME <- minimize(ttOUTCOME, include = "?", details=TRUE, show.cases=TRUE, dir.exp = "AWILL, AACTION, COMPLEX")
isOUTCOME

# produce table
stargazerSol(results = isOUTCOME, outcome = "INFL",
             type = "html", show.cases=TRUE, out = "isOUTCOME.html")



# XY plot for solution


pimplot (data =  mydata,
          outcome = "INFL" ,
          results = isOUTCOME,
          all_labels = TRUE,
         jitter = TRUE,
         fontsize = 6)


# identify simplifying assumptions 

SAINFL <- psOUTCOME$SA
SAINFL

# export SAs as table
write.csv(SAINFL, "SAINFL.csv")

# easy counterfactuals

ECINFL <- isOUTCOME$i.sol$C1P1$EC
ECINFL

# export ECs as table
write.csv(ECINFL, "ECINFL.csv")

#identify typical and deviant cases for SMMR

#typical
typINFL <- smmr(results = isOUTCOME,
                outcome = "INFL" ,
                match = FALSE ,
                cases = 1,
                max_pairs = 5)
typINFL


# deviant coverage

devcovINFL <- smmr(results = isOUTCOME,
                outcome = "INFL" ,
                match = FALSE ,
                cases = 4,
                term = 1,
                max_pairs = 7)
devcovINFL

#pair deviant coverage--IIR case for comparison (identify omitted conjunction)

pairINFL <- smmr(results = isOUTCOME,
                   outcome = "INFL" ,
                   match = TRUE ,
                   cases = 4,
                   term = 1,
                   max_pairs = 2)
pairINFL



####### ~INFL (REPORTED IN PAPER)########
ttoutcome <- truthTable(data=mydata, outcome = "~INFL", 	
                        conditions = "AWILL, AACTION, COMPLEX, CONTEST", 
                        incl.cut=0.75, sort.by="incl, n", complete=TRUE, show.cases=TRUE)
ttoutcome

# truth table for INFL, with dcc only
ttoutcome1 <- truthTable(data=mydata, outcome = "~INFL", 	
                         conditions = "AWILL, AACTION, COMPLEX, CONTEST", 
                         incl.cut=0.75, sort.by="incl, n", complete=TRUE, show.cases=TRUE, dcc=TRUE)

ttoutcome1 

# The first row has three cases of which 1 dcc. equally row 1 (the third row) has 1 dcc case out of three cases, but consistency is low.

# we exclude the third row (nr. 1) because of low consistency

# We set the raw consistency threshold at .8, to include the first two rows 
ttoutcome <- truthTable(data=mydata, outcome = "~INFL", 	
                        conditions = "AWILL, AACTION, COMPLEX, CONTEST", 
                        incl.cut=0.8, sort.by="incl, n", complete=TRUE, show.cases=TRUE)
ttoutcome

# export truth table

write.csv(ttoutcome$tt, "ttinfl.csv")

#conservative solution

csoutcome <- minimize(ttoutcome, details=TRUE, show.cases=TRUE, use.tilde=TRUE)
csoutcome

# produce table

stargazerSol(results = csoutcome, outcome = "~INFL",
             type = "html", show.cases=TRUE, out = "csoutcome.html")

# parsimonious solution
psoutcome <- minimize(ttoutcome, include="?", details=TRUE, show.cases=TRUE )
psoutcome

# produce table

stargazerSol(results = psoutcome, outcome = "~INFL",
             type = "html", show.cases=TRUE, out = "psoutcome.html")

# intermediate solution, Enhanced Standard Analysis. We expect that ceteris paribus, ~AWILL, ~AACTION,  ~COMPLEX would lead to ~INFLUENCE. There is no expectation for CONTEST.
# exclude untenable assumptions and perform ESA

#identify contradictory simplifying assumptions

CSA <- LR.intersect(psOUTCOME, psoutcome)
CSA

# there are no necessary conditions to be considered and no impossible remainders
# truth table for enhanced standard analysis

ettoutcome <- esa (ttoutcome ,
                   contrad_rows = c(CSA))
ettoutcome

# enhanced intermediate solution
eisoutcome <- minimize(ettoutcome, include = "?", details=TRUE, show.cases=TRUE, dir.exp = "~AWILL, ~AACTION, ~COMPLEX")
eisoutcome

# produce table

stargazerSol(results = eisoutcome, outcome = "~INFL",
             type = "html", show.cases=TRUE, out = "eisoutcome.html")


# XY plot for solution


pimplot (data =  mydata,
         outcome = "~INFL" ,
         results = eisoutcome,
         all_labels = TRUE,
         jitter = TRUE,
         fontsize = 6)

# identify simplifying assumptions 

SAinfl <- psoutcome$SA
SAinfl

# export SAs as table
write.csv(SAinfl, "SAnotinfl.csv")

# easy counterfactuals

ECinfl <- eisoutcome$i.sol$C1P1$EC
ECinfl

# export ECs as table
write.csv(ECINFL, "ECnotinfl.csv")

#identify typical and deviant cases for SMMR

# typical

typinfl <- smmr(results = eisoutcome,
                outcome = "~INFL" ,
                match = FALSE ,
                cases = 1,
                max_pairs = 5)
typinfl

# deviant coverage

devcovinfl <- smmr(results = eisoutcome,
                   outcome = "~INFL" ,
                   match = FALSE ,
                   cases = 4,
                   term = 1,
                   max_pairs = 7)
devcovinfl

# deviant consistency
devconinfl <- smmr(results = eisoutcome,
                   outcome = "~INFL" ,
                   match = FALSE ,
                   cases = 3,
                   term = 1,
                   max_pairs = 7)
devconinfl


# typical-deviant consistency

pairsconinfl <- smmr(results = eisoutcome,
                     outcome = "~INFL" ,
                     match = TRUE ,
                     cases = 3,
                     term = 1,
                     max_pairs = 7)
pairsconinfl

#pair deviant coverage--IIR case for comparison (identify omitted conjunction)

pairinfl <- smmr(results = eisoutcome,
                 outcome = "~INFL" ,
                 match = TRUE ,
                 cases = 4,
                 term = 1,
                 max_pairs = 2)
pairinfl

######ROBUSTNESS TEST 1: LESS RESTRICTIVE ANALYSIS FOR ~INFL--NOT ADOPTED IN PAPER######

# truth table for ~INFL, with all cases
rttoutcome <- truthTable(data=mydata, outcome = "~INFL", 	
                        conditions = "AWILL, AACTION, COMPLEX, CONTEST", 
                        incl.cut=0.75, sort.by="incl, n", complete=TRUE, show.cases=TRUE)
rttoutcome

# truth table for INFL, with dcc only
rttoutcome1 <- truthTable(data=mydata, outcome = "~INFL", 	
                         conditions = "AWILL, AACTION, COMPLEX, CONTEST", 
                         incl.cut=0.75, sort.by="incl, n", complete=TRUE, show.cases=TRUE, dcc=TRUE)

rttoutcome1 

# The first row has three cases of which 1 dcc. equally row 1 (the third row) has 1 dcc case out of three cases, but consistency is low.

# we include the third row (nr. 1) because UNESCO_ALE is a unique case of INFL, with extremely low salience of the issue. Also, PRI is nevertheless above 0.5 and thus this row is more sufficient for ~INFL than for INFL

# We set the raw consistency threshold at .668, to include the first three rows 
rttoutcome <- truthTable(data=mydata, outcome = "~INFL", 	
                        conditions = "AWILL, AACTION, COMPLEX, CONTEST", 
                        incl.cut=0.668, sort.by="incl, n", complete=TRUE, show.cases=TRUE)
rttoutcome

# export truth table

write.csv(rttoutcome$tt, "ttinfl.csv")

#conservative solution

rcsoutcome <- minimize(rttoutcome, details=TRUE, show.cases=TRUE, use.tilde=TRUE)
rcsoutcome

# produce table

stargazerSol(results = rcsoutcome, outcome = "~INFL",
             type = "html", show.cases=TRUE, out = "csoutcome.html")

# parsimonious solution
rpsoutcome <- minimize(rttoutcome, include="?", details=TRUE, show.cases=TRUE )
rpsoutcome

# produce table

stargazerSol(results = rpsoutcome, outcome = "~INFL",
             type = "html", show.cases=TRUE, out = "psoutcome.html")

# intermediate solution, Enhanced Standard Analysis. We expect that ceteris paribus, ~AWILL, ~AACTION,  *COMPLEX would lead to ~INFLUENCE. There is no expectation for CONTEST.
# exclude untenable assumptions and perform ESA

#identify contradictory simplifying assumptions

rCSA <- LR.intersect(psOUTCOME, rpsoutcome)
rCSA

# there are no necessary conditions to be considered and no impossible remainders
# truth table for enhanced standard analysis

rettoutcome <- esa (rttoutcome ,
                   contrad_rows = c(rCSA))
rettoutcome

# enhanced intermediate solution
reisoutcome <- minimize(rettoutcome, include = "?", details=TRUE, show.cases=TRUE, dir.exp = "~AWILL, ~AACTION, ~COMPLEX")
reisoutcome
isOUTCOME
# produce table

stargazerSol(results = reisoutcome, outcome = "~INFL",
             type = "html", show.cases=TRUE, out = "reisoutcome.html")


# XY plot for solution


pimplot (data =  mydata,
         outcome = "~INFL" ,
         results = reisoutcome,
         all_labels = TRUE,
         jitter = TRUE,
         fontsize = 6)

# identify simplifying assumptions 

rSAinfl <- rpsoutcome$SA
rSAinfl

# export SAs as table
write.csv(rSAinfl, "rSAnotinfl.csv")

# easy counterfactuals

rECinfl <- reisoutcome$i.sol$C1P1$EC
rECinfl

# export ECs as table
write.csv(rECINFL, "rECnotinfl.csv")

##this less restrictive variant is not adopted in the paper because a) it produces more model ambiguities than the restrictive model,and  b) it comes with more inconsistencies.

#####ROBUSTNESS ANALYSIS#######

# due to the small-N character of our analysis we do not assess changes of the frequency cutoff. Moreover, since our sets were calibrated directly from the qualitative information (i.e. there is no raw dataset) and we have high conceptual certainty about these decisions, robustness against calibration decisions is also not assessed.

####changes in raw consistency thresholds####

###for INFL#####

####sensitivity range
# create object containing the conditions
conds <- c("AWILL", "AACTION", "COMPLEX", "CONTEST")

# change the threshold by value 0.01 and do so 40 times


rob.inclrange(data=mydata, step=0.01, max.runs=40, outcome = "INFL", conditions=conds, incl.cut = 0.830, n.cut = 1, include = "?" )

ttOUTCOME

# The initial solution for INFL is somewhat very sensitive to changes: the Boolean formula stays the same if the raw consistency threshold is placed between 0.51 and 0.83 (there are no actual truth table rows within this range). In other words, we set the threshold as low as is possible while maximizing coverage. This makes us confident that we took the right decision. At the same time, and unsurprisingly, the formula would change if we set a higher RC threshold.


#### fit-oriented robustness

# initial solution
IS <- isOUTCOME

# we test against 2 alternative raw consistency thresholds. Both of these alternatives only make limited conceptual sense though: with fuzzy sets, it is not necessary to restrict RC to 1 (TS1). Moreover, a threshold of 0.5 is normally considered much too low (TS2), and rightly so as in the respective tt row, 2 out of 3 cases are contradictory.The robustness test, in short, should simply be seen as indicative and taken with caution.We are confident in our choice of the raw consistency threshold, as it is the one that covers most typical cases while minimizing the existence of DCCK cases.

# TS1: require consistency of 1
TS1ttOUTCOME <- truthTable(data=mydata, outcome = "INFL", 	
                               conditions = "AWILL, AACTION, COMPLEX, CONTEST", 
                               incl.cut=1, sort.by="incl, n", complete=TRUE, show.cases=TRUE)
TS1ttOUTCOME

TS1 <- minimize(TS1ttOUTCOME, include = "?", details=TRUE, show.cases=TRUE, dir.exp = "AWILL, AACTION, COMPLEX")
TS1
  
# TS2: include row 3 with 1 more case of INFL and RC of 0.5

TS2ttOUTCOME <- truthTable(data=mydata, outcome = "INFL", 	
                           conditions = "AWILL, AACTION, COMPLEX, CONTEST", 
                           incl.cut=0.5, sort.by="incl, n", complete=TRUE, show.cases=TRUE)
TS2ttOUTCOME

TS2 <- minimize(TS2ttOUTCOME, include = "?", details=TRUE, show.cases=TRUE, dir.exp = "AWILL, AACTION, COMPLEX")
TS2


# create test set
TS <- list(TS1, TS2)

# calculate fit-oriented robustness

RF <- rob.fit(test_sol = TS, initial_sol=IS, outcome="INFL")
RF

# the intermediate solution for INFL is not very robust regarding parameters its coverage--that is, by changing the raw consistency threshold, we further increase / decrease the coverage. However, it is reasonably, though not perfectly robust regarding consistency. The solution reported in the paper coincides little with the minimal test set, but a reasonably well with the maximal test set. 

#### case-oriented robustness
rob.xyplot(test_sol = TS, initial_sol = IS, outcome = "INFL", all_labels = TRUE, fontsize = 3.5, jitter = TRUE, area_lab = TRUE)

# the robustness rank of 4 indicates again low robustness: cases could change all possible kinds of status when different raw consistency thresholds are chosen. In the upper left quadrant, we have three possible cases that would be newly covered by the alternative solutions. On the lower right, we see five shaky cases since their status as typical or deviant would change with these changes. Conversely, 7 cases on the diagonal are part of the robust core that would not change. only 40 per cent of the typical cases are robust. 

# obtain names of cases & Boolean expressions
rob.cases(test_sol = TS, initial_sol = IS, outcome = "INFL")

###for ~INFL####

IS <-eisoutcome

####sensitivity ranges

# change the threshold by value 0.01 and do so 40 times


rob.inclrange(data=mydata, step=0.01, max.runs=40, outcome = "~INFL", conditions=conds, incl.cut = 0.8, n.cut = 1, include = "?" )

ttoutcome

# The initial solution for ~INFL is somewhat very sensitive to changes: the Boolean formula stays the same if the raw consistency threshold is placed between 0.67 and 0.8 (No tt rows there). At the same time, and unsurprisingly, the formula would change if we set a higher RC threshold.


#### fit-oriented robustness####

# initial solution
IS <- eisoutcome

# determining test set
ettoutcome

# we test against 1 alternative raw consistency threshold--the same as under robustness test 1, to include 2 more cases of INFL.

# TS1: raw consistency of 0.668
TS1ttoutcome <- truthTable(data=mydata, outcome = "~INFL", 	
                           conditions = "AWILL, AACTION, COMPLEX, CONTEST", 
                           incl.cut=0.668, sort.by="incl, n", complete=TRUE, show.cases=TRUE)
TS1ttoutcome

TS1 <- minimize(TS1ttoutcome, include = "?", details=TRUE, show.cases=TRUE, dir.exp = "~AWILL, ~AACTION, ~COMPLEX")
TS1


# create test set
TS <- list(TS1)

# calculate fit-oriented robustness

RF <- rob.fit(test_sol = TS, initial_sol=IS, outcome="~INFL")
RF

# the intermediate solution for ~INFL is  very robust regarding parameters of fit:  it perfectly overlaps with the robust core. At the same time, the set coincidence between the initial solution and the minimal and maximal test set is quite low.This means that when changing the raw consistency thresholds, other cases would be included that the initial solution does not cover. 

# case-oriented robustness
rob.xyplot(test_sol = TS, initial_sol = IS, outcome = "~INFL", all_labels = TRUE, fontsize = 3.5, jitter = TRUE, area_lab = TRUE)

# the robustness rank of 2 indicates again reasonable robustness: there are no shaky cases, but there are possible cases, i.e. cases that would be additionally covered with the test solution. These are two cases in the upper left. Other than that, most cases cluster on the diagonal, indicating high robustness. 

# obtain names of cases & Boolean expressions
rob.cases(test_sol = TS, initial_sol = IS, outcome = "~INFL")
